The information presented here is placed in the public domain, and was written by Doug Burke. The notebook used to create this page is available, and questions can be asked using the GitHub issues page or via Twitter: https://twitter.com/doug_burke.

It will probably make a lot-more sense if you have already read the previous entry in this poorly-defined series (if you can call two notebooks a series).

What am I up to now?

The whole point of the first notebook was to try out the cool-looking units package, motivated by the paper "Experience Report: Type-checking Polymorphic Units for Astrophysics Research in Haskell" by Takayuki Muranushi and Richard A. Eisenberg (Haskell Symposium 2014, Gothenburg, Sweden). However, I couldn't install the package on the system-installed version of ghc on my laptop, so in the interest of people who may be trying to follow along at home, I stuck with version 7.6.3 of ghc and used the dimensional-tf package instead. As I do have newer versions of ghc on this laptop, I can also try out units at the possible expense of making it harder for people to replicate the notebook:

  • switch to version 7.8.4 of ghc

  • in a new directory, create a new cabal sandbox; e.g.

    % cabal sandbox init
    % cabal install ihaskell units-defs integration --dry-run
    % cabal install ihaskell units-defs integration
    ...
  • remove the existing ihaskell configuration directories

    % rm -rf ~/.ihaskell
    % rm -rf ~/.ipython/profile_haskell
    
    

    I did this because without it the IHaskell kernel would fail with the message

    IHaskell: Bad interface file: /home/dburke/local/ghc-7.8.4/lib/ghc-7.8.4/base-4.7.0.2/Prelude.hi
        mismatched interface file versions (wanted "7063", got "7084")
  • start the IHaskell notebook with

    % mkdir notebooks
    % cabal exec IHaskell -- notebook --serve=notebooks/

Last time the notebook was run


In [1]:
import Data.Time
getCurrentTime


2015-03-04 00:56:25.676914 UTC

Load the modules

I need to turn on a few GHC extensions to use the units code:


In [2]:
-- TypeFamilies is needed when defining Unit instances
:set -XTypeFamilies

-- TypeOperators is needed for defining the EnergyFlux type
:set -XTypeOperators

-- DataKinds is needed for defining AngleDim and Radians
:set -XDataKinds

Unlike last time, I'm not going to create any plots, so I don't need the Chart modules.

I do want the integration routine, and I also have to load in the relevant modules from units and units-defs, which provides the definitions of the SI units I need.


In [3]:
-- Numerics
import Numeric.Integration.TanhSinh

In [4]:
-- Units
import Data.Metrology
import Data.Metrology.SI
import Data.Metrology.Show

The units package provides two "forms" of the API: a general version (the Data.Metrology.Poly/Data.Metrology.SI.Poly modules), which provides more control, and a "simpler" version, which I use here. The difference is not in the function names, but in their types, with the polymorphic version being more general than the simpler, monomorphic, version. The code below would be unchanged if I used the .Poly versions of the modules, I just would have to supply more complex types (or at least, that's what I believe; I haven't actually tried it ;-).

Using the units module

As with dimensional-tf, there's a way to convert a number into a type with units. In this case it is %, but the values to the right are now Haskell types, rather than values as they are in dimensional-tf, which means that they are handled in a different manner. The units package provides syntax for "combining" units; here :/ is used to indicate division, but there's also :* for multiplication, :^ for raising to a power, and :@ for combining terms (there's examples of most, if not all, of these later on in the notebook).

I start by defining the speed of light, which can be done by sepcifying the numeric value and then the units, separated by the % operator. This is similar to the *~ operator of dimensional-tf, except that the units are defined by types rather than values (Haskell starts type names with a capital letter). As we shall see, there's not much difference between these two approaches for the code I write here, but that's not really a fair comparison as you don't get to see all the mis-steps I made along the way, and this is not a complex piece of code that needs much maintenance or extension.


In [5]:
c = 299792458.0 % (Meter :/ Second)
:type c


c :: forall n. Fractional n => Qu '['F Length One, 'F Time ('P 'Zero)] 'DefaultLCSU n

The inferred type for this value is as equally baroque (for the non-Haskeller) as the dimensional-tf version, which for reference was

c = 299792458.0 P.*~ (P.meter P./ P.second)
:type c
c :: forall a. Fractional a => Quantity (Dim (S Z) Z (N (S Z)) Z Z Z Z) a

In this case the first argument to Qu, which is the term within `[ ... ], indicates the dimensions - so here length of 1 and time of -1 (Peano numbers again being used to indicate values in the type, with aliases for common numbers, such as One), with the second term being a mysterious DefaultLCSU value (which I am going to ignore, since I can get away with writing this notebook without having to use it, as it only really comes into play when using the polymorphic version), and the third is the type of the numeric value.

As with dimensional-tf, there are more-readable types which can be specified. In this case Velocity, which is defined as Length :/ Time:


In [6]:
c :: Velocity
c = 299792458.0 % (Meter :/ Second)

:type c


c :: Velocity

While % promotes a number into a unit, # extracts a value from a value in the specified units. These are similar in spirit to the *~ and /~ operators of dimensional-tf, but they have rather eye-glazing types:


In [7]:
:type (%)
:type (#)


(%) :: forall n unit. (Subset (CanonicalUnitsOfFactors (UnitFactorsOf unit)) (CanonicalUnitsOfFactors (LookupList (DimFactorsOf (DimOfUnit unit)) 'DefaultLCSU)), Subset (CanonicalUnitsOfFactors (LookupList (DimFactorsOf (DimOfUnit unit)) 'DefaultLCSU)) (CanonicalUnitsOfFactors (UnitFactorsOf unit)), Unit unit, UnitFactor (LookupList (DimFactorsOf (DimOfUnit unit)) 'DefaultLCSU), Fractional n) => n -> unit -> Qu (DimFactorsOf (DimOfUnit unit)) 'DefaultLCSU n
(#) :: forall n unit. (Subset (CanonicalUnitsOfFactors (UnitFactorsOf unit)) (CanonicalUnitsOfFactors (LookupList (DimFactorsOf (DimOfUnit unit)) 'DefaultLCSU)), Subset (CanonicalUnitsOfFactors (LookupList (DimFactorsOf (DimOfUnit unit)) 'DefaultLCSU)) (CanonicalUnitsOfFactors (UnitFactorsOf unit)), Unit unit, UnitFactor (LookupList (DimFactorsOf (DimOfUnit unit)) 'DefaultLCSU), Fractional n) => Qu (DimFactorsOf (DimOfUnit unit)) 'DefaultLCSU n -> unit -> n

To combine values, you have to use the operators provided by the module; these use a naming scheme where '|' appears to the side that a quantity is. This means that to multiply c by 0.99 you could say either of


In [8]:
0.99 *| c
c |* 0.99


2.9679453342e8 m/s
2.9679453342e8 m/s

To multiply two values together you use |*|, which also keeps track of the combined units:


In [9]:
c |*| c


8.987551787368176e16 m^2/s^2

or add together values with consistent units, which also shows how compound units such as cm or kg are formed using the :@ operator to combine the types:


In [10]:
(4 % Meter) |+| (3 % (Centi :@ Meter))


4.03 m

I used brackets to make it clearer what terms were related, but fortunately the precedence of the various operators means that you can often get away without using them. In this particular case no brackets are needed since:


In [11]:
4 % Meter |+| 3 % Centi :@ Meter


4.03 m

As you'd hope, you get errors when trying to combine values with incompatible units, as shown in the next set of examples:


In [12]:
4 % Meter |+| 0.03


Couldn't match type ‘'F Data.Dimensions.SI.Length One : Normalize '[]’ with ‘'[]’
In the expression: 4 % Meter |+| 0.03
In an equation for ‘it’: it = 4 % Meter |+| 0.03

In [13]:
0.99 * c


Couldn't match type ‘'F Data.Dimensions.SI.Length One : Normalize '['F Data.Dimensions.SI.Time ('P 'Zero)]’ with ‘'[]’
In the expression: 0.99 * c
In an equation for ‘it’: it = 0.99 * c

In [14]:
1 % Meter |+| 1 % Second


Couldn't match type ‘'F Data.Dimensions.SI.Length One : Normalize '['F Data.Dimensions.SI.Time ('P 'Zero)]’ with ‘'[]’
In the expression: 1 % Meter |+| 1 % Second
In an equation for ‘it’: it = 1 % Meter |+| 1 % Second

If you use the normal mathematical operators then you also get errors, even if the units match:


In [15]:
(1 % Meter) + (1 % Meter)


Couldn't match type ‘'[]’ with ‘'['F Data.Dimensions.SI.Length One]’
In the first argument of ‘(+)’, namely ‘(1 % Meter)’
In the expression: (1 % Meter) + (1 % Meter)
In an equation for ‘it’: it = (1 % Meter) + (1 % Meter)

Couldn't match type ‘'[]’ with ‘'['F Data.Dimensions.SI.Length One]’
In the second argument of ‘(+)’, namely ‘(1 % Meter)’
In the expression: (1 % Meter) + (1 % Meter)
In an equation for ‘it’: it = (1 % Meter) + (1 % Meter)

Putting this together, I can calculate the value of $\frac{1}{2} m v^2$ for $m = 1~{\rm kg}$ and $v=c$:


In [16]:
answer = 0.5 *| (1 % Kilo :@ Gram) |*| c |*| c
:type answer
answer


answer :: Qu '['F Mass One, 'F Length ('S ('S 'Zero)), 'F Time ('P ('P 'Zero))] 'DefaultLCSU Double
4.493775893684088e16 (kg * m^2)/s^2

This can be given an explicit type, rather than letting the compiler infer one:


In [17]:
answer2 :: Energy
answer2 = answer
answer2


4.493775893684088e16 (kg * m^2)/s^2

Using an explicit type catches dimensional errors in your formula, but you may have to use the redim function - as shown below - to ensure that the types can be compared: redim has a compile-time cost, when it forces the units to be compared, but there is no run-time cost to using it. It is needed because the type system can't easily match up units where only the order of the dimensions differs.


In [18]:
-- This is going to fail because I have missed out multiplying by a velocity,
-- so the units of the calculation are just MASS * LENGTH / TIME, which
-- does not match energy. If redim is not included then you get, for this
-- case, a more-verbose error, but there is the possibility in this case that
-- it's just a limitation of the type system rather than an actual units
-- mis-match.
answer3 :: Energy
answer3 = redim (0.5 *| (1 % Gram) |*| c)


Couldn't match type ‘'F Data.Dimensions.SI.Length ('P 'Zero) : Normalize '['F Data.Dimensions.SI.Time ('S 'Zero)]’ with ‘'[]’
In the expression: redim (0.5 *| (1 % Gram) |*| c)
In an equation for ‘answer3’: answer3 = redim (0.5 *| (1 % Gram) |*| c)

As previoulsy mentioned, # is used to extract the value in the given units (so it is similar in spirit to the /~ operator of dimensional-tf). For this case it works whether an explicit type was given (answer2) or not (answer):


In [19]:
answer # Joule
answer2 # Mega :@ Joule


4.493775893684088e16
4.493775893684088e10

The current version of units (well, really, it's units-defs) does not contain some of the types I want to use (e.g. Parsec and Astronomical Unit). I can use the existing functionality to allow me to create a scale value such as


In [20]:
parsec :: Length
parsec = 3.085678e16 % Meter

which I can then use in equations


In [21]:
3 *| parsec


9.257034e16 m

In [22]:
(3 *| parsec) # Meter


9.257034e16

However, it's not a type, so I can not convert a length into parsecs using # (well, not without manually converting the answer using the correct scaling factor), as shown below (the error message is so large because I have made a logical mistake here, using the wrong "thing" on the right-hand side of the # operator, and the type system is therefore really confused):


In [23]:
(12.3e18 % Meter) # parsec


Couldn't match type ‘DimFactorsOf (DimOfUnit (Qu '['F Data.Dimensions.SI.Length One] 'DefaultLCSU Double))’ with ‘'['F Data.Dimensions.SI.Length One]’
In the expression: (12.3e18 % Meter) # parsec
In an equation for ‘it’: it = (12.3e18 % Meter) # parsec

Could not deduce (Subset (CanonicalUnitsOfFactors (UnitFactorsOf (Qu '['F Data.Dimensions.SI.Length One] 'DefaultLCSU Double))) '[Meter],
Elem Meter (CanonicalUnitsOfFactors (UnitFactorsOf (Qu '['F Data.Dimensions.SI.Length One] 'DefaultLCSU Double))))
arising from a use of ‘#’
from the context (Fractional n) bound by the inferred type of it :: Fractional n => n at :1:1-26
In the expression: (12.3e18 % Meter) # parsec
In an equation for ‘it’: it = (12.3e18 % Meter) # parsec

Could not deduce (Unit (Qu '['F Data.Dimensions.SI.Length One] 'DefaultLCSU Double)) arising from a use of ‘#’
from the context (Fractional n) bound by the inferred type of it :: Fractional n => n at :1:1-26
In the expression: (12.3e18 % Meter) # parsec
In an equation for ‘it’: it = (12.3e18 % Meter) # parsec

It is possible, however, to create a type for parsecs. Following the units tutorial I came up with the following, which relates Parsec to Meter.


In [24]:
data Parsec = Parsec

instance Unit Parsec where
  type BaseUnit Parsec = Meter
  
  -- Note the _ argument!
  conversionRatio _ = 3.085678e16

The conversionRatio function has the following type:


In [25]:
:type conversionRatio


conversionRatio :: forall unit. Unit unit => unit -> Rational

which says that the first argument must be an instance of the Unit type class (which is what the instance Unit ... section above does), and returns a Rational value (which is a particular representation of a floating-point number). The documentation for units makes it clear that the first argument to this function must never be used (technically, I believe it's that this value should never be evaluated); to do this we use the "placeholder" syntax of _, which tells Haskell to ignore the argument. For users not used to Haskell, it might seem strange to ignore an argument to a function, or to not "evaluate" it. The reason for ignoring this argument is that it's actually the type of the argument that is important, rather than its value. The type tells the compiler which particular version of conversionRatio to use (it's a bit like how interfaces in Python or Java are used to select what code to use, but please don't hold strongly to this viewpoint as thinking about objects is likely to keep you confused when reading Haskell), and so it doesn't need to use the actual value. In such situations, it is common for the value to be set to undefined, which is a special value, in that

  • I don't really want to try and explain it here

  • it can be set to any type (i.e. it's polymorphic)

  • evaluating the value is going to cause an error (hence the requirement that user code should never, ever, be used)

Here's an example of such an error, which isn't too informative, and you won't be seeing elsewhere in this notebook:


In [26]:
undefined


Prelude.undefined

With all this I can now create a length in parsecs:


In [27]:
bob = (3::Double) % Parsec

The type, shown below, indicates that it's a length. I could have said

bob :: Length
bob = (3::Double) % Parsec

which provides more-readable code, and error messages, but is the same as the version without the explicit type. The restriction of 3 to a Double type is just to simplify the types a little bit (it removes the forall n. Fractional a => constraint on a type, since Double is an instance of the Fractional type class).


In [28]:
:type bob
:type (3::Double) % Meter


bob :: Qu '['F Length One] 'DefaultLCSU Double
(3::Double) % Meter :: Qu '['F Length One] 'DefaultLCSU Double

Its value is reported in metres, since that is the base unit of the Parsec type:


In [29]:
bob


9.257034e16 m

With this type, I can also convert other lengths into parsecs. The two forms below are the same, since mega can be considered to be a short form for saying Mega :@, and the () are not needed, since the % operator takes precedence over #, but I have left them in to make it clearer what is going on.


In [30]:
(2.2e28 % Meter) # mega Parsec


712971.3469778765

In [31]:
(2.2e28 % Meter) # Mega :@ Parsec


712971.3469778765

I can do the same for the Hubble Constant, that is define a type, but it doesn't feel as sensible to me, since we do not really talk about units of km/s/Mpc in quite the same way as we do parsecs.


In [32]:
-- Hubble constant in units of km/s/Mpc
data HubbleConstant = HubbleConstant

instance Unit HubbleConstant where
  type BaseUnit HubbleConstant = Hertz
  
  -- As shown in later examples, this could have used something
  -- like (1 % (Kilo @: Meter :/ Mega @: Parsec :/ Second))
  -- to calculate the conversion ratio, but I wanted to
  -- use an explicit value here:
  conversionRatio _ = 1e3 / (1e6 * 3.085678e16)

With this, I can easily create a value equal to 70 km/sMpc:


In [33]:
h70 = 70 % HubbleConstant
h70
:type h70


2.268545194929607e-18 s^-1
h70 :: forall n. Fractional n => Qu '['F Time ('P 'Zero)] 'DefaultLCSU n

Alternatively, a value can be created just by specifying the units directly, with the same result (i.e. the type is the same):


In [34]:
h70 = 70 % (Kilo :@ Meter :/ Second :/ Mega :@ Parsec)
h70
:type h70


2.268545194929607e-18 s^-1
h70 :: forall n. Fractional n => Qu '['F Time ('P 'Zero)] 'DefaultLCSU n

One advantage to having the HubbleConstant type is that I can convert a frequence into km/s/Mpc, as shown below, but I personally don't feel this is particularly compelling:


In [35]:
(7.2e-19 % Hertz) # HubbleConstant


22.2168816

For this notebook, an alternative to creating the HubbleConstant unit would be to create a function which adds in the correct units, such as the hubbleConstant function below:


In [36]:
hubbleConstant x = x % (Kilo :@ Meter :/ Second :/ Mega :@ Parsec)
hubbleConstant 70
:t hubbleConstant


2.268545194929607e-18 s^-1
hubbleConstant :: forall n. Fractional n => n -> Qu '['F Time ('P 'Zero)] 'DefaultLCSU n

As mentioned before, it seems to be personal preference for which to prefer (at least for the tasks I perform in this notebook), as the results are the same:


In [37]:
hubbleConstant 72.4
72.4 % HubbleConstant


2.3463238873271936e-18 s^-1
2.3463238873271936e-18 s^-1

The angular-diameter and luminosity distances

It is now possible to calculate the angular-diameter and luminosity distances. As with the dimensional-tf versions, I have included explicit types to ensure that the inputs and outputs are as expected, since without the signatures the compiler would infer very-generic types. A complication is in having to remember which variant of the operators should be used (e.g. *| or |*), but the compiler complains quite loudly if you get it wrong (another reason for having the signatures). As previously mentioned, redim is essentially doing book-keeping at compile time; I didn't seem to need to include them here, but as they have no run-time cost it seems a sensible thing to include.


In [38]:
-- angular-diameter and luminosity distances (units of the Hubble length)

dcH om z = result (absolute 1.0e-6 (trap f 0 z))
  where
    -- unlike the previous notebook, I have made f a local definition
    -- here, so it does not need to be sent in the om value as
    -- it is already in scope
    f z = let ol = 1 - om
              t = (1 + z)^2 * (1 + om*z) - (2 + z) * ol * z
          in 1 / sqrt t
         
daH om z = dcH om z / (1 + z)
           
dlH om z = dcH om z * (1 + z)

-- now the versions which return values in physical units

da :: Double -> Frequency -> Double -> Length 
da om h0 z = redim (daH om z *| c |/| h0)

dl :: Double -> Frequency -> Double -> Length
dl om h0 z = redim (dlH om z *| c |/| h0)

Now it's time to compare the results to those from dimensional-tf (I do not expect any significant numerical differences, what I am interested in is how to write and use the values). I use the hubbleConstant function here, but the result would have been the same if I had used 70 % HubbleConstant instead.


In [39]:
da 0.3 (hubbleConstant 70) 2 # Mega :@ Parsec


1726.6206914696234

The next step is to calculate the transverse distances, by multiplying the angular-diameter distance da by an angle in radians. The dimensional-tf package has in-built support for angles, but units does not. For this section I treat angles as simple numbers (so do not take advantage of any dimensional support), but do add in a simple type to represent angles, to better document the intent of the routines. I try creating a Radian unit at the end of this notebook to compare to this approach.

Here I create an Angle type which stores a single number as a Double, which I am going to assume is in radians. To help enforce this, I also create several helper routines which will create an Angle type given several different input values (e.g. a value in degrees or arc seconds). Since this is being done in a single IHaskell notebook it is quite easy to subvert the safety provided by these routines by creating an Angle value manually. In "real" Haskell code, the Angle data type would be defined in a separate module which would only allow you to create values using the routines such as radian, rather than exposing the ability to create an Angle type directly.


In [40]:
-- The Angle type has a single numeric value, which is taken to be the
-- value in radians.

data Angle = Angle { toRadians :: Double }
  deriving Show
  
-- Convert a number in radians, degrees, arcminutes, arcseconds to an Angle data type

radian :: Double -> Angle
radian x = Angle x

degree :: Double -> Angle
degree x = radian (x * pi / 180)

arcmin :: Double -> Angle
arcmin x = degree (x / 60)

arcsec :: Double -> Angle
arcsec x = arcmin (x / 60)

As a check, here's the output for these routines:


In [41]:
radian 1
degree 1
arcmin 1
arcsec 1


Angle {toRadians = 1.0}
Angle {toRadians = 1.7453292519943295e-2}
Angle {toRadians = 2.908882086657216e-4}
Angle {toRadians = 4.84813681109536e-6}

The value from an Angle can be retrieved using the toRadians function:


In [42]:
:type toRadians
toRadians (degree 1)


toRadians :: Angle -> Double
1.7453292519943295e-2

With these, the adist routine can be written as:


In [43]:
adist :: Double -> Frequency -> Double -> Angle -> Length
adist om h0 z angle = da om h0 z |* toRadians angle

Since the fourth argument has a type of Angle, I can not call it with a number, otherwise the compiler complains:


In [44]:
adist 0.3 (hubbleConstant 70) 2 1


No instance for (Num Angle) arising from the literal ‘1’
In the fourth argument of ‘adist’, namely ‘1’
In the expression: adist 0.3 (hubbleConstant 70) 2 1
In an equation for ‘it’: it = adist 0.3 (hubbleConstant 70) 2 1

Supplying an actual Angle value allows me to compare to the dimensional-tf version; as above, the numeric values are essentially the same


In [45]:
adist 0.3 (hubbleConstant 70) 2 (arcsec 1) # Kilo :@ Parsec


8.370893333112804

In [46]:
adist 0.3 (hubbleConstant 70) 2 (arcmin 2.4) # Mega :@ Parsec


1.205408639968244

And now for the luminosity distance:


In [47]:
dl 0.3 (hubbleConstant 70) 2


4.7950159338113435e26 m

In [48]:
dl 0.3 (hubbleConstant 70) 2 # Mega :@ Parsec


15539.58622322661

Perhaps surprisingly, units-defs does not come with a definition of a year, which means that calculating this distance in units of Giga light years is going to require some set up. Although it is not clear what a good general-purpose definition is, I am going to define one following the IAU recommendation, which also happens to be the version that dimension-tf used.


In [49]:
-- Following the dimensional package, define a year as 365.25 days,
-- which is an IAU recommendation 
-- http://www.iau.org/science/publications/proceedings_rules/units/
-- and who am I to argue with the demoters of Pluto?
--
data Year = Year

-- I could have made the base unit be second and use
-- a conversion ratio of 31557600, but I wanted to see
-- how this sort of "chaining" of units worked.
instance Unit Year where
  type BaseUnit Year = Hour
  conversionRatio _ = 365.25 * 24

-- Light years
data LightYear = LightYear

instance Unit LightYear where
  type BaseUnit LightYear = Meter
  
  -- How efficient is the following: is the number calculated at compile time
  -- or evaluated at each conversion? This approach is in contrast to the
  -- Unit instance of HubbleConstant, where I explicitly included
  -- the numeric values.
  --
  -- Without the realToFrac function, the answer would be a Fractional,
  -- but conversionRatio returns a Rational, so I have to do an
  -- explicit type cast here.
  --
  conversionRatio _ = realToFrac (c # Meter :/ Year)

After all that, I can finally calculate the distance in light years:


In [50]:
dl 0.3 (hubbleConstant 70) 2 # Giga :@ LightYear


50.68335841199911

Calculating luminosities

The next step is to convert from flux to luminosity. This requires a non-SI unit, erg, which is provided by Data.Units.CGS:


In [51]:
import Data.Units.CGS (Erg(..))

With this I can create energies and convert between Ergs and Joules:


In [52]:
1 % Erg

1 % Joule |/| 1 % Erg


1.0e-7 (kg * m^2)/s^2
1.0e7

To make the type signature cleaner I invent the following type to represent fluxes:


In [53]:
type EnergyFlux = Power %/ Area

The luminosity is therefore just (again, ignoring the K correction term):


In [54]:
luminosity :: Double -> Frequency -> Double -> EnergyFlux -> Power
luminosity om h0 z flux = redim (flux |*| area)
  where
    r = dl om h0 z
    area = 4 * pi *| r |*| r

Unfortunately I can't seem to create fluxes on the fly, since this causes type errors because of the order of the units, as shown in the following error message (which turns out to be because of a missing call to redim for the flux rather than a problem with the units):


In [55]:
luminosity 0.3 (hubbleConstant 70) 2 (2e-13 % (Erg :/ Centi :@ Meter :^ sTwo :/ Second))


Couldn't match type ‘Data.Dimensions.SI.Length’ with ‘Data.Dimensions.SI.Time’
In the fourth argument of ‘luminosity’, namely ‘(2e-13 % (Erg :/ Centi :@ Meter :^ sTwo :/ Second))’
In the expression: luminosity 0.3 (hubbleConstant 70) 2 (2e-13 % (Erg :/ Centi :@ Meter :^ sTwo :/ Second))
In an equation for ‘it’: it = luminosity 0.3 (hubbleConstant 70) 2 (2e-13 % (Erg :/ Centi :@ Meter :^ sTwo :/ Second))

One way around this is to define a helper function which includes a call to redim (or to include the redim call directly in the original call):


In [56]:
ergFlux :: Double -> EnergyFlux
ergFlux x = redim (x % (Erg :/ Second :/ Centi :@ Meter :^ sTwo))

With this definition, I can create a luminosity


In [57]:
lum = luminosity 0.3 (hubbleConstant 70) 2 (ergFlux 2e-13)
lum


5.778564550704302e38 (kg * m^2)/s^3

and convert it to the units I require:


In [58]:
lum # Watt
lum # (Erg :/ Second)


5.778564550704302e38
5.778564550704302e45

The critical density

The final part of the original notebook was to calculate the critical density of the Universe as a function of Hubble's constant, using the equation:

$$\rho_{\rm crit} = \frac{3 H^2}{8 \pi G}$$

where $G = 6.67428 \times 10^{-11}\ {\rm m}^3 / \rm{kg} / \rm{s}^2$:


In [59]:
criticalDensity :: Frequency -> Density
criticalDensity h = redim (3 *| h |*| h |/| (8 * pi *| bigG))
  where
    -- bigG is only defined local to the criticalDensity routine
    bigG = 6.67428e-11 % (Meter :^ sThree :/ Kilo :@ Gram :/ Second :^ sTwo)

This lets me calculate the same value as with dimensional-tf:


In [60]:
c70 = criticalDensity (70 % HubbleConstant)
c70


9.203899006458817e-27 kg/m^3

To finish off I need to create a unit for Astronomical Units:


In [61]:
-- From http://en.wikipedia.org/wiki/Astronomical_unit 
data AstronomicalUnit = AstronomicalUnit

instance Unit AstronomicalUnit where
  type BaseUnit AstronomicalUnit = Meter
  conversionRatio _ = 149597870700

With this, I can convert to the (completely made up) units of kilo-tonnes per cubic au:


In [62]:
c70 # (Kilo :@ Ton :/ (AstronomicalUnit :^ sThree))


30.814000174159233

Representing angles

The following is what I came up with when attempting to add in support for an "angle dimension". I make no guarantee that it's sensible, but it seems to work.


In [63]:
-- This is based on the Number and Dimensionless types from Data.Metrology.Units 
data AngleDim = AngleDim

instance Dimension AngleDim where
  type DimFactorsOf AngleDim = '[]

type instance DefaultUnitOfDim AngleDim = Radian

data Radian = Radian
instance Unit Radian where
  type BaseUnit Radian = Canonical
  type DimOfUnit Radian = AngleDim
  type UnitFactorsOf Radian = '[]

data Degree = Degree
instance Unit Degree where
  type BaseUnit Degree = Radian
  conversionRatio _ = realToFrac (pi / 180.0)

data ArcMinute = ArcMinute
instance Unit ArcMinute where
  type BaseUnit ArcMinute = Degree
  conversionRatio _ = (1/60)

data ArcSecond = ArcSecond
instance Unit ArcSecond where
  type BaseUnit ArcSecond = ArcMinute
  conversionRatio _ = (1/60)

-- I have already defined an Angle, so I add Type to the name
-- here to have a different symbol name.
type AngleType = MkQu_U Radian

With this, I can create a value in radians and convert it to different units:


In [64]:
angle = 0.4 % Radian
angle


0.4

As this is a dimensionless number, the first argument to Qu is empty:


In [65]:
:type angle


angle :: forall n. Fractional n => Qu '[] 'DefaultLCSU n

In [66]:
3 *| angle
:t 3 *| angle


1.2000000000000002
3 *| angle :: forall n. Fractional n => Qu '[] 'DefaultLCSU n

Since the angles are dimensionless, you do not need to use the special operators on them, as shown below. I shall keep using the operators from units just to act as a form of documentation.


In [67]:
3 * angle
:t 3 * angle


1.2000000000000002
3 * angle :: forall n. Fractional n => Qu '[] 'DefaultLCSU n

Values can be converted between units:


In [68]:
(3 *| angle) # Degree
(1.2 % ArcMinute) # ArcSecond


68.7549354156988
72.0

The units package defines a dimensionless type, Number, which can be used to "extract" the numeric value from these angles.


In [69]:
(3 *| angle) # Number
:t (3 *| angle) # Number


1.2000000000000002
(3 *| angle) # Number :: forall n. Fractional n => n

With all this I can write a version of adist that uses such an angle:


In [70]:
adist2 :: Double -> Frequency -> Double -> AngleType -> Length
adist2 om h0 z angle = da om h0 z |*| angle

The results are the same as above:


In [71]:
adist2 0.3 (hubbleConstant 70) 2 (1 % ArcSecond) # Kilo :@ Parsec


8.370893333112804

In [72]:
adist2 0.3 (hubbleConstant 70) 2 (2.4 % ArcMinute) # Mega :@ Parsec


1.2054086399682438

I did find one "surprising" thing with this version: even though the fourth argument is given as AngleType, you can actually give a number, which the compiler will interpret as a value in radians. So, repeating the calculation for one arcsecond, I get the same result, after explicitly converting the angle to degrees):


In [73]:
adist2 0.3 (hubbleConstant 70) 2 (pi/3600/180) # Kilo :@ Parsec


8.370893333112804

Did I learn anything?

Overall, I would say that there's not a huge difference in using the two packages dimensional-tf and units (I best hurry and add a disclaimer pointing out that this really isn't enough code to really determine a preference, and also that I don't want this to be an "a is better than b" style post).

The two packages do things differently, which has implications on how things are defined, what needs to be written, how readable are the error messages, and what is done in the value system as compared to the type system, but for the purposes of this notebook the actual code implementing the equations are not too dissimilar:

  • dimensional-tf
bigG = 6.67428e-11 P.*~ (P.metre P.^ P.pos3 P./ P.kilo P.gram P./ P.second P.^ P.pos2)
criticalDensity h = P._3 P.* h P.^ P.pos2 P./ (P._8 P.* P.pi P.* bigG)
  • units
bigG = 6.67428e-11 % (Meter :^ sThree :/ Kilo :@ Gram :/ Second :^ sTwo)
criticalDensity :: Frequency -> Density
criticalDensity h = redim (3 *| h |*| h |/| (8 * pi *| bigG))

As dimensional-tf is the older package, it has many of the units I needed, and it was easy to add new ones, such as parsec. The units and units-defs packages are newer, and do not contain as many values or types, but it is also not hard to add in new ones, as shown above.

The units package avoids the need for the qualified terms that I used with the dimensional-tf version, for instance the P.XXX forms above, but this is in part due to defining its own mathematical operators, such as |* and |/|. The units version is slightly-more readable to me, but this is not enough code to really decide if either has a strong advantage. There's also the fact that you benefit most from the things that don't get shown here - i.e. all the wrong turns and errors that the type system complains about until you get it right - and I really have not had enough experience of the two packages to make any comment on that point.

One concern, when using a units package, is that of numeric precision: if the values all get converted to some base unit for each dimension, could large (or small) values lead to numerical issues. Since I used Double types here, and the values didn't get too huge (in part because I never got around to calculating cosmological volumes), this didn't seem to be a problem. Section 5.2 of "Experience Report: Type-checking Polymorphic Units for Astrophysics Research in Haskell" describes how units can handle this; I'm not sure how to cope with this in the dimensional-tf package.

The end

There you go; I hope you enjoyed it. If you have any questions, then please use the GitHub issues page or contact me on Twitter at https://twitter.com/doug_burke.